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This paper introduces the method of anisotropic diffusion oscillation reduction (ADOR) for shock 
wave computations. The connection is made between digital image processing, in particular, image 
edge detection, and numerical shock capturing. Indeed, numerical shock capturing can be formulated 
on the lines of iterative digital edge detection. Various anisotropic diffusion and super diffusion oper- 
ators originated from image edge detection are proposed for the treatment of hyperbolic conservation 
laws and near-hyperbolic hydrodynamic equations of change. The similarity between anisotropic 
diffusion and artificial viscosity is discussed. Physical origins and mathematical properties of the 
artificial viscosity is analyzed from the kinetic theory point of view. A form of pressure tensor is 
derived from the first principles of the quantum mechanics. Quantum kinetic theory is utilized to 
arrive at macroscopic transport equations from the microscopic theory. Macroscopic symmetry is 
used to simplify pressure tensor expressions. The latter provides a basis for the design of artificial 
viscosity. The ADOR approach is validated by using (inviscid) Burgers' equation in one and two 
spatial dimensions, the incompressible Navier-Stokes equation and the Euler equation. A discrete 
singular convolution (DSC) algorithm is utilized for the spatial discretization. 







I. INTRODUCTION 



Shock wave is a common phenomenon in nature, such as in aerodynamics and hydrodynamics. MathematicaUy, 
nonhnear hyperbohc conservation equations provide a good description to shock waves. The construction of numerical 
schemes that are capable of shock capturing for hyperbolic and inviscid hydrodynamic equations is a major objective 
of computational methodology. However, it is noted that the concept of discontinuity does not apply in digital 
computations. Therefore, a shock in computational sense may refer to systematic, rapid variation of function values 
over a few grid points. The difficulty is that hyperbolic equation may have weak solutions that are discontinuous at 
the so-called shock front. Such a discontinuity will cause Gibbs' oscillations in a high (spatial) order numerical scheme. 
The numerically induced oscillations are usually amplified in (time) iterations. A variety of numerical schemes have 
been proposed for shock capturing. As early as 50 years ago, a solution to this problem was constructed by von 
Neumann and Richtmyer |l| . The essence of their approach is to introduce small artificial viscosity so that a smooth 
solution can alway be attained in a finite difference approximation. A variety of modifications to von Neumann and 
Richtmyer's method are made in the last 5 decades to address problems of possible failure in spatial scaling and errors 
due to additional momentum flux and production, as well as unbalanced heat flux and its over production in the 
simulation of hydrodynamic conservation laws. 

A different approach was proposed by Godunov Q to construct a full solution by using low order piecewise discon- 
tinuous approximations. Such a piecewise solution is a good approximation at the smooth regions, and is capable of 
representing the shock front over a small region of grid. The knowledge of wave propagation and wave interaction 
is built in the numerical scheme in the form of a Riemann solver. The Godunov's approach has been extended to 
higher-order schemes [^j-^. In doing so, a higher-order approximation is constructed in the smooth region, while near 
the shock the solution is still of first order accuracy. The Godunov method is very stable and thus, easy to design 
and use However the major disadvantage of this method is the complexity introduced into a numerical scheme 
through a Riemann solver. Such a drawback reduces its computational efficiency. 

Another general approach is the hybrid scheme, which utilizes a high-order scheme for a smooth region while 
using a low-order scheme near a discontinuity. A linear combination of these two types of schemes is then used 
at each interface using weight factors which may be nonlinear functions of the local flow field. DeBar constructed 
a linear hybridization of first-order and second-order difference schemes as early as 1968 0. Harten and Zwas |^ 
devised another early linear hybridization scheme. Boris and Book reported a blending algorithm which yields sharp 
discontinuities without oscillations 1^. A total variation diminishing (TVD) scheme was proposed by Harten [Tc| ] 
to control the spurious oscillations in the numerical solution by using the total variation as a measure. The TVD 
scheme typically degenerates to first-order accuracy at locations with smooth extrema and was later generalized to 
an essentially non-oscillatory (ENO) scheme p]-^. The major idea of the ENO scheme is to suppress spurious 
oscillations near the shock or discontinuities, while maintaining a higher-order accuracy at smooth regions. This line 
of thinking was further polished recently in a weighted essentially non-oscillatory (WENO) scheme [Q. The WENO 
approach takes a linear combination of a number of high-order schemes of both central difference and up-wind type. 
The central difference type schemes have a larger weight factor at the smooth region while the up-wind schemes play 
a major role at the shock or discontinuity. In general, these approaches are quite expensive since checks are performed 
before making a decision at each grid point. In terms of accuracy, all existing methods are at best of first order near 
the shock or discontinuity. 

Perhaps modified artificial viscosity methods are the most popular approaches in practical computations. Unfor- 
tunately, the artificial viscosity smears shocks over three or more grid zones, which can lead to serious errors in the 
physical interpretation of the numerical results. Special care is required to ensure that the smeared numerical shock is 
consistent with the true thickness of the shock in a practical problem under study. This difficulty has led to enormous 
and continuous effort at developing efficient and robust approaches. One approach is to locally refine the computa- 
tional grid 1^ . An alternative approach is the use of an adaptive mesh . Both methods are aimed at matching 
between the physical shock and the numerically smoothed one with respect to the spatial extension. However, a major 
drawback is their restriction for extremely small time step sizes as required by the Courant constraint. In many cases, 
they have to be formulated in an implicit scheme, which imposes an extra complexity in practical implementation 
and an extra requirement in computer memory. Another problem is the unbalanced momentum flux and production, 
and additional heat flux and heat production in hydrodynamic equations. 

The addition of a viscous term to the inviscid hydrodynamic equations or the hyperbolic equation has a physical 
justification, that the true physical flow has no discontinuities. Therefore, mathematical model can be modified so 
as to refiect the true physics. In their original work, von Neumann and Richtmyer ^ have introduced an artificial 
viscosity qn r of the form 

QNR - (V • v)' (1) 



1 



for the equation of motion. Their term is of second order in gradient of the velocity field v (The velocity gradient 
is a second rank tensor and has the divergence of velocity as its component) and will not be very sensitive to small 
gradients. Landshoff proposed an additional term that vanishes less rapidly for small gradients 



A generic form that contains both qnr and (i.e., aV • v + /3 (V-v) ) has been carefully studied by many researchers 



Jl9|-p2[. In particular, Caramana et al discussed a number of intuitive criteria for this form and the parameteri- 
zation of artificial viscosity. Special considerations are given to dissipativity, Galilean invariance, self-similar motion 
invariance, wave front invariance and viscous force continuity. This generic form, as it was proposed for hydrodynamic 
conservation laws, has its advantages over other approaches for certain fluid mechanical computations. However, it 
is very unstable for certain hyperbolic equations due to the Courant-Friedrich-Lewy (CFL) stability condition. For 
example, it does not work well in resolving a sharp shock front for inviscid Burgers' equation with a smooth initial 
condition. We believe that this failure is due to the high-order in the gradient of the Neumann and Richtmyer form, 
which is too sensitive to large gradients at a sharp shock front. Despite the fact that a number of functional forms 
of the artificial viscosity have been proposed, the procedure is still ad hoc. For a given form, parameter selection is 
quite tricky and often varies from problem to problem. 

The existence of so many different approaches for shock capturing indicates both the importance and the difficulty 
of the problem. The objective of the present paper is to introduce an anisotropic diffusion oscillation reduction 
(ADOR) approach for shock wave computations. The method of anisotropic diffusion in association with a partial 
differential equation was proposed by Perona and Malik ||2^ for digital image processing in 1990. Since then, much 
research has been stimulated in image processing and applied mathematics communities ||2^-|3^. The method uses 
the heat equation, which has a solution of the Gaussian type, as a low pass filter to eliminate noise in an image, 
while it detects and preserves the edges of the image. It has been shown that the Perona-Malik equation provides a 
computational approach to image segmentation, noise removal, edge detection, and image enhancement. In a recent 
work, a generalized Perona-Malik equation was proposed for image restoration and edge enhancement |^^. In fact, 
digital edge detection has much in common with numerical shock capturing. In this work, we introduce a set of 
generalized anisotropic diffusion operators and edge enhancing functionals for shock wave computations. 

This paper is organized as follows: Section II is devoted to kinetic theory analysis of the equation of change and 
artificial viscosity. The physical origin and mathematical properties of artificial viscosity are discussed from the 
kinetic theory point of view. A set of hydrodynamic equations describing fiuid flow consisting of microscopic particles 
is derived from the quantum Boltzmann equation, i.e., the Waldmann- Snider equation [ ^J37| . The latter can be 
regarded as a consequence of a reduction of the von Neumann equation, or the quantum Liouville equation to the 
level of a single particle density operator. The mathematical form of the pressure tensor is analyzed by using the 
group theory consideration. A set of generalized artificial viscosities, artificial heat and artificial virial correction are 
proposed as the results of the kinetic theory analysis. Generalized Perona-Malik equation is discussed in Section III 
for shock wave computations. This approach follows a very different line of thinking as it was originally proposed for 
image analysis and processing. An edge-controlled image enhancing functional is introduced for shock representation. 
Anisotropic diffusion and diffusion super operators are introduced for fast and effective shock capturing. Numerical 
experiments are presented in Section IV. A few standard test problems, including Burgers' equation and inviscid 
Burgers' equation in one- and two-dimensions, the incompressible Euler and Navier-Stokes equations are used for 
exploring usefulness, testing efficiency and illustrating the validity of the ADOR approach. A recently developed 
discrete singular convolution (DSC) algorithm p8|-p[ is utilized for the numerical integration. This paper ends with 
a conclusion. 



Artificial viscosity was original proposed for shock capturing in association with hydrodynamic equations in the 
fiuid mechanics. The choice of its form and parameter should be consistent with and motivated by the physical origin 
of hydrodynamic equations. Since quantum theory is the foundation for the modern fiuid mechanics, a microscopic 
analysis is presented for appropriate understanding from the kinetic theory. The latter, such as the Boltzmann 
equation, serves as a theoretical basis for hydrodynamic conservation laws. 



V-v. 



(2) 




II. KINETIC THEORY ANALYSIS OF ARTIFICIAL VISCOSITY 
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A. Microscopic Analysis 



Consider the flow of a quantum gas system consisting of total N particles in a volume S. Its behavior is governed 
by the Schrodinger equation 



at 



(3) 



where i?'^^ is the self-adjoint Hamiltonian of the system and ^' is a vector in the Hilbert space associated with 
the system. For the description of physical observable, we adopt the density operator p^^^ whose time evolution is 
governed by the quantum Liouville equation 



' dt P ^ h 



(4) 



A physical observable O^^^ is a Hermitian operator of the Hilbert space and has the expectation value given by 

<0W >=Tri,.,^OWpW, (5) 
where Tri^...^7v is the trace over all the states of the N particle system and in particular, 



< iW Tn... 



= 1 



(6) 



gives the normalization of the total probability for finding all of the N particles in the volume S. Here 1*^^^ is the 
identity operator of the system. A standard form for the Hamiltonian is given by 



(7) 



where Ki is the kinetic energy (including internal state energy) operator of particle i and Vij is the potential operator 
of particles i and j. The physical behavior of the A'^-particle system is far too complicated to compute by any means 
as a macroscopic gas flow may consist of 10^'^ particles or more. Fortunately, for ordinary gases, it is sufScient to 
consider the state of a typical particle, say particle 1 



Pi 



«=iVTr2,..,^pW. 



In general, the state of n particles is defined 



PV^,„ = N{N - 1) ■ • • (iV - n + l)Tr„+i,...,jvP 



(JV) 



The time evolution of particle 1 is governed by 



(1) 



L'- 'Pi + ir2, 1/12 P12 ^ [ '^1 



dt 



(1) 



1 r 
h 



Vl2,p?^ 



(8) 



(9) 



(10) 



where H^^^ = Ki. This is the first member of the BBGKY hierarchy [ ^ , ^ in the quantum form. The general form 
of the BBGKY hierarchy E|,|l is given by 



■ "Pl,-,n _ f^{n)M ^ ^,{n+l) (ri+1) 



(11) 



where V[^.t'^2^i is the potential superoperator between particles 1, • • • , n and particle n + 1. This set of equations 

is formal and exact, since no approximation has been made. However, to determine the time evolution of p\^\ it is 

(2) . . . 

necessary to know the behavior of , which is, in turn, determined by the second order BBGKY equation and the 
latter involves three-particle density operator ^223- 
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B. Mesoscopic analysis 



Kinetic theory attempts to explain macroscopic properties in terms of microscopic properties of the atoms and/or 
molecules based on the classical or quantum mechanics. Typical macroscopic observations deal with iV(~ 10^^) 
particles over a volume much larger than the size of the individual molecules, and over a time period much longer 
compared to the time scale of the individual molecular dynamics. An ab-initio description of the time dependence 
of such macroscopic observations from the quantum Liouville equation is practically impossible not only because of 
the conceptual difficulty with the meaning of measurement, but also due to the large number of degrees of freedom 
involved. The kinetic theory approach simplifies the iV-particle problem dramatically by looking at the behavior of 
one typical particle under the influence of all other particles. 

The most successful kinetic theory has been based on the Boltzmann equation ji^js^ , ^ , ^ . Specifically this has 
been related to the hydrodynamic equations and has given rise to molecular expressions for the transport coefficients. 
Various attempts have been made to understand the Boltzmann equation from the A^-body Liouville equation, or 
utilizing the BBGKY hierarchy described in the last subsection. A first principle "derivation" of the Boltzmann 



equation from the Liouville equation has been achieved by Bogoliubov | |43| , and independently by Green 46 1. The 
Wang Chang-Uhlenbeck-de Boer equation |^7|] was the first attempt to generalize the quantum Boltzmann equation, 
and to account for the internal degrees of freedom of a polyatomic gas. It was not generally realized until 1957 that 
polyatomic gases can only be treated properly by a quantum mechanical method due to the fact that the internal 
energy levels of such molecules are, in general, degenerate. Waldmann p6[ | and, independently. Snider [ |37[ | introduced 
a quantum kinetic equation, the Waldmann-Snider equation, 

_ (1) (1) (2)(^ (1) (1) 

where ^Li2 the M0ller superoperator for quantum mechanical scattering 

Ql^, - lim e<^^e-<^ K (13) 

t — ^—00 

The Waldmann-Snider equation is still the basis for the kinetic theory of quantum gas fiow [EspS]. Note that a 
dramatically simplified version of the Boltzmann equation provides the (lattice) Boltzmann gas (LBG) approach for 
the computational fiuid dynamics in the recent years. 



C. Macroscopic analysis 

The most important properties of a fluid flow are physical observables, or physical measurements. For conservative 
observables, the equations of change of physical observables give rise to conservation laws. The equations of change 
are important not only because they govern the time dependencies of the macroscopic quantities but also because 
they are needed for solving kinetic equations using the Chapman-Enskog method ||5^] . The equations of change for 
physical observables are derived in this subsection. 

A single-particle observable for particle 1 is denoted by tpi. The expectation value of the physical observable ^pl is 

determined by the singlet density operator p'"^^ according to 

(^i) = TriV'ipi''. (14) 

The most important physical observables are the number density ni, stream velocity vi and kinetic energy . Firstly, 
the physical observable for the number density ni is the Dirac delta distribution V'l — Si = S{r — Vi) so that 

ni = {Si). (15) 

Secondly, the stream velocity vi is given by 

vi = — Tri-^(pi5i + Sipi)pi, (16) 
ni Zmi 

where mi is the mass of particle 1. Finally, the kinetic energy ef per particle is given by 

ef = — Tn (J—{pl5i + 2pi • Sipi + SipD - ^rmvlSi] pi. (17) 
ni \8toi 2 J 
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These expressions are appropriately symmetrized to account for the fact that position and momentum do not commute, 
while a physical observable is a Hermitian operator. Although this operator approach is used here for for the number 
density ni, stream velocity Vi and kinetic energy ef^, phase space expressions can be easily obtained for these fluid 
dynamic quantities by using the standard Wigner representation. The equations of change are of course independent 
of the detailed method used for the expression of various quantities. Equations of change for a particle observable is 
obtained from the quantum Boltzmann equation according to 

d{ipi)i .dipi 1 ^) 1 (2) (1) (1) , , 



In case that a physical observable tpi is time independent, the second term on the left hand side of Eq. ( |18D vanishes. 
Conservation laws are the equation of continuity 



the equation of motion 



and the kinetic energy equation 



^ = -V-(niVi); (19) 



9vi 

ninii-^ + nimivi • Vvi = —V • Pi, (20) 



anief __ _^ _ ^^^^^^^ ^ ^ ^^^^^^ _ p, ^^^^ ^ ^21) 



dt 

where the superscript t denotes the transpose and the kinetic heat flux is 



rj. \ I Pi \ (pi -mivi)% 
qi = Tri <^ vi di } pi, (22) 



^mi / 2mi 

where { }s designates appropriately operator-symmetrized quantities. Here the collisional heat flux is 



1 



qcou = -z Tri2 / di^ <S 



8 m 



-1 



I"! + r2 ly 

^^12 \ „ (1) (1) 



X (pi + p2 - 2toiVi) . \ ^Il.A 'pI ' (23) 



and the kinetic energy production results from the collisional transfer of energy 

a^ = -^Tri2|(5i+^2)(pi-p2^ 



8mi 1^ ari 

(Pi-P2)l ^L.Jl^p'^i''- (24) 



9ri 

The pressure tensor Pf^ has kinetic and collisional contributions 

Pi = Pf + Pl°^^ (25) 

where the kinetic pressure tensor is 

pf = 7^('^iPiPi + Pi^iPi + (Pi<5iPi)* + PiPi<5i)i - nimivivi (26) 
4toi 

and the collisional pressure tensor is given by 

Pi°" = -^Tri2(ri-r2)Viyi2 

' 5 i ^ + ^ (n " r2) - v\ d.n,J^^p^i\ (27) 



1 



2 2 
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The transpose P* of the pressure tensor Pi enters to couple the convective energy to the internal kinetic energy per 
particle , heat flux contributions arise from kinetic qf^ and coUisional qcoii motion and finally there is a production 
term which involves the transfer between kinetic and potential energies. Although the equations of continuity and 
motion are consistent with the conservation of the number of particles and total momentum in the absence of bound 
pairs, the presence of the production term in the kinetic energy implies that the total kinetic energy is not conserved 
because of the possible conversion to potential energy due to the non-locality of the collisions (in the macroscopic 
sense). Thus it is necessary to look at the equation of change for the potential energy per particle eY 

meX = ^Tri2((5i + 52)Vi2p''il (28) 
This is obtained in a manner consistent with pair particle interactions and shown to be of the form 

^ = -V.(niVier + qr)-^^, (29) 

where 

= 7^Tri2 [(pi - miVi)Ji + 5i(pi - miVi)] Vi2nL,J^\^2^ ■ (30) 
4m 1 

The production of potential and kinetic energy exactly cancels and there is an added heat flux contribution c\( 
associated with the conductive potential energy flow. 



D. Pressure tensor and artificial viscosity 



The estimation of transport coefficients from the Boltzmann equation is quite complicated and involves many 
aspects. The most important transport quantities for physical and engineering applications are the mass, momentum 
and energy. Sometimes angular momentum transport may also be important for certain physical phenomena and 
thus an equation of change for angular momentum may be required. However, angular momentum conservation is 
rarely discussed in the context of hydrodynamic conservation laws. This is because the angular momentum is not 
measured like other physical observables in experiments. It is noted that all hydrodynamic conservation laws have the 
structure of convection, conduction and production. The equation of momentum conservation leads to shock wave as 
the divergence of the pressure tensor vanishes. The existence of shock wave devastates the numerical simulation as 
noted by non Neumann and Richtmyer Q. They introduced artificial viscosity of the form of (V • v)^ to overcome 
the numerical difficulty. In fact, a much better choice can be selected by analyzing the form of the pressure tensor. 

The Pressure tensor is a tensor of second rank and has two indices. It has contributions from the kinetic transport 
P^ and collision transport P'^"" (Here, the subscript 1 is omitted for convenience). For example, the kinetic transport 
P^y is the rate of transport (kg m/s)/(m^s) in the x-direction, of the momentum with respect to the y-direction. The 
hydrodynamic pressure 

Peq = nksTU, (31) 

given by the equation of state, is the local equilibrium contribution to the pressure tensor. Here, U is the identity 
tensor of second rank and fcs is the Boltzmann constant. In general, the gas is not at local equilibrium and the 
structure of the nonequilibrium part of the pressure tensor, 

n = P - Peq (32) 

can be analyzed according to the irreducible representation of the improper three-dimensional rotational group 0(3). 
As such, n can be expressed as a linear combination the irreducible representations 

n = HU + e-P'' + n(2\ (33) 

where scalar, antisymmetric and symmetric traceless parts of the second rank tensor are given by 

H4u:n 
P'' ^ \e U 

n(2) = l[n + n*] - ^u[U:n]. (34) 
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Here e is the Levi-Civita tensor, H is a scalar (one dimension), P° is a vector (three dimensions) and II^^^ is a tensor 
(five dimensions). In the theory of the Boltzmann equation and irreversible thermodynamics, the pressure tensor is 
treated as proportional to velocity gradient as in the case of a classical Newtonian flow. The velocity gradient, Vv, 
is also a second rank tensor and can also be decomposed into three different components of order zero, one and two, 
similar to Eq. (||), 

Vv = iu:Vv 

V X V = -^e-Vv 

[Vv](') =i[Vv + (Vv)*]-iuV.v. (35) 

Here, Vv is also a nine-dimensional quantity. In principle, there are 81 possible coefficients connecting the pressure 
tensor and the velocity gradient. Symmetry analysis indicates that there are three phenomenological equations relating 
the corresponding components of 11 and Vv 

n = -»y°v-v 

= -vy X V 

n(2) =_2r,0[Vv](2), (36) 

where riy,rfj. and rf are the viscosities of bulk, rotational and shear. 

The pressure tensor given by Eq ( p6[ ) provides an excellent approximation for a gas system near the local equilibrium, 
i.e. the Boltzmann distribution 

For the purpose of numerical stability and shock wave simulation, we can impose artificial pressure tensor as 

Hart = -Cl^V • V 
P:rt = -C"V X V 

ni?) = -2C°[Vv](2), (38) 

where C^y , and C^^ are artificial bulk viscosity, artificial rotational viscosity and artificial shear viscosity, respectively. 
These artificial viscosities modify the non equilibrium part of the pressure tensor. In fact, from the point of view of 
computations, the equilibrium part of the pressure tensor can also be modified, i.e. a term which is proportional to 

P„T = enfcsTU. (39) 

PnT can be interpreted either as artificial heat or artificial virial correction. Hence, we have the artificial pressure 
tensor of the form 

Part = PnT + Ha,, U + e-P^,, + U'^^l (40) 

= enksTU - Cv^ ■ vU - Cr^-V x v - 2C"[Vv](2). (41) 

It should be noted that there is an advantage in keeping only a part of the artificial pressure tensor that corresponds 
to the non-equilibrium pressure tensor. The latter makes no additional contribution to conservative quantities, such 
as, number density, momentum and energy, as required by the Frcdholm alternative for the existence of a solution for 
the kinetic equation in the Chapman-Enskog method |Q . 

However, for a gaseous system far from the local equilibrium, higher-order terms in the velocity gradient become 
important. Large velocity gradient can build up from special boundary condition, such as the boundary layer phenom- 
ena, or from the lack of relaxation in a very dilute gas flow. Therefore artificial viscosities and artificial heat of Eqs. 
( ^ ) can be modified to include higher-order velocity gradient contributions. Recognizing that the pressure tensor is 
a second rank tensor, thus all velocity gradients contribute in one of the forms of H, P'^ and II'^^ In this regard, the 
artificial viscosity form proposed by von Neumann and Richtmyer [Q, (V • v)^, contributes to H. Certainly, artificial 
viscosities of the forms of V • v, V x v and [VvJ^^^^ can be used for numerical computations. All of the aforementioned 
forms are at least of the first order in gradient, which may lead a strong smoothing at the edge of a shock front and 
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thus lead to large errors in numerical applications. Furthermore, these forms can result in instability when the shock 
front is very sharp, such as in inviscid Burgers' equation. For these reasons, it is appropriate to propose a general 
expression for the artificial viscosity and artificial heat so as to allow all coefficients to be functions of V • v, V x v 
and [Vv](2) 



e 



V.v,||Vxv||,||[Vv](2)||) (42) 
C{^-Cl^fv-v,||Vxv||,||[Vv](2)||) (43) 



C° = C°(v-v,||Vxv||,||[Vv](^)||) (44) 
C°-C"fv-v,||Vxv||,||[Vv]W||) (45) 



where || • || denotes the magnitude and is computed as 



l|A|| = ^5:Af (46) 

for a vector A, and 



m = JT.^l (47) 

for a tensor B. Obviously, artificial viscosity of von Neumann and Richtmyer (V • v)^, becomes a special case of 
the present treatment. The von Neumann and Richtmyer form can be classified either as an artificial heat term, PnT, 
or as an artificial bulk viscosity term, IlartU. 

The inclusion of an artificial rotational viscosity can be an efficient way for the treatment of fast rotational fiow. 
It is noted that the rotational viscosity is in the original Navier-Stokes equation derived from the kinetic theory. 
However, for flows with natural boundary conditions, the angular momentum is also conserved and it does not couple 
the linear momentum except through internal spinor relaxations. However, this is no longer the case in a flow with 
an irregular geometry. For such a case, the internal-conversion between angular momentum and linear momentum is 
substantial and therefore the rotational viscosity, should be retained in both theoretical modeling and numerical 
simulation. A detailed discussion and numerical simulation of this aspect will be accounted elsewhere. In the next 
section, we analyze the shock capturing algorithm further from the point of view of image processing, particularly, 
image edge detection. 

III. OSCILLATION REDUCTION BY ANISOTROPIC DIFFUSION 

Although anisotropic diffusion was original proposed for image processing, there is much in common between 
digital image processing and computational fluid dynamics. An image function /(r) is a two-dimensional projection of 
certain physical quantities, such as matter, velocity, energy, electromagnetic field, etc., under appropriate illumination 
conditions. The edges in an image usually refer to rapid changes in some physical properties, such as geometry, 
illumination, and reflectance. Mathematically, a discontinuity may be involved in the function representing such 
physical properties. Therefore image edges are very similar to shocks in fluid dynamics. Numerical shock capturing 
can be formulated on the lines of iterative digital edge detection. Edge detection is a key issue in image processing, 
computer vision, and pattern recognition. A variety of algorithms, such as the Sobel operator, the Prewitt operator, 
the Canny operator |5l| and the DSC algorithm are proposed for image edge detection and representation. 
Anisotropic diffusion is a promising new mathematical algorithm for image edge detection and image processing. This 
basic idea can be adopted and modified for numerical shock capturing in association with the hyperbolic conservation 
laws. The Perona-Malik algorithm ||2^] is reviewed before the method of anisotropic diffusion oscillation reduction 
(ADOR) is discussed. 

A. The Perona-Malik equation 

The basic idea behind the Perona-Malik algorithm is to evolve an original image, /(r), under an edge-controlled 
diffusion operator pj 
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^^=V.[d(||Vu(r,i)||)V^(r,<)] 
"(r,0)=/(r). (48) 

Here, (i(||Vu||) is a generalized diffusion coefficient which is so designed that its values are very small at the edges of 
an image. Many edge stopping functions (i(| Vw|) are appropriate for anisotropic diffusion. For example, the Gaussian 

d(||Vu||) = e-l^"l'/2"' (49) 

and the Lorcntz 

d(||V^.||) - _} ^ (50) 

are both suitable for edge representation. They provide perceptually similar results in practical applications. Numer- 
ically, the diffusion coefficient becomes very small near an image edge due to the effect of edge stopping functions 
d(l|Vw||). As a result, the image edge is preserved in the diffusion process. The pixel values at a non-edge part will be 
smoothed and reduced due to the substantial diffusion coefficient prescribed by (i([|Vw||). Perona and Malik argued 
that the solution of their anisotropic diffusion equation has no additional maxima (minima) which does not belong to 
the initial image data. However, this point has been challenged recently |^,^. It is well-known that this anisotropic 
diffusion algorithm may break down when the gradient generated by noise is comparable to image edges and features. 
Numerically, this can be alleviated by using a regularization procedure. 



B. Shock capturing by anisotropic diffusion 

Hyperbolic conservation laws of the type 

^ + V.F(.)^0 (51) 

describe the rate of change of a physical quantity u given by the generalized convection V • F(7i). Without the 
balance of conduction and/or production, Eq. (|5l|) may have a discontinuous solution. The task is to construct a 
stable computational scheme which is capable of resolving the "shock" . We first note that computationally, if a shock 
is defined as a discontinuity, there is no shock to capture. This is because, the original notion of discontinuity is 
undefined on a discrete mesh. Therefore, a numerical shock is characterized by rapid variation of function values over 
a small grid zone. Numerical shock capturing, at best, is globally a first order approximation scheme for resolving 
the large gradient feature of the true solution. However, locally, it is preferred to compute the solution as accurate as 
possible, except at the discontinuity of the true solution. To this end, we introduce an anisotropic diffusion term to 
Eq. (|1|) 

+ V • Fiu) = V • [d,{\\\7u{r, t)\\)Vu{r,t)] (52) 

where the diffusivity, (ii(||VM(r, t)||, is chosen such that it is essentially zero except at a numerical shock position. 
Obviously, ||Vu(r, i)|| is important for shock detection. Apparently, Eq. ( |5^ ) reduces to the Perona-Malik equation 
( ^ ) without the convection term. However, the selections of the anisotropic diffusivity in these two equations are 
entirely different and serve opposite purposes. For example, one may choose (ii(||Vii(r, i)|| as 

di = d}ln[(Vii)2 + 1], (53) 

where d\ is a constant. Equation (p3) differs much from the Gaussian and the Lorentz. 



The derivation of Eq. 



has its roots in the conservation of a physical quantity involving a phenomcnological 



fiux and its divergence. However, as a computational algorithm, an anisotropic diffusion of the form 

V ■ F{u) = ri(||V?/(r,t)||)V2u(r, t) (54) 



dt 

can be used. In fact, expressions in both Eqs. ( |5^ ) and ( p^ are efficient for numerical shock capturing. Here, Fi is 
chosen to smooth the oscillations near a shock and is essentially zero in other regions. 
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C. Edge enhancing functional and super diffusion 



From the hydrodynamic point of view, the governing equation of a conservative quantity has a general structure, 
i.e., the rate of change is balanced by convection, conduction and production. Therefore, the nonlinear advective 
motion can be counterbalanced by an appropriate production. Therefore, we propose a real-valued, bounded shock 
(edge) enhancing functional 

e(||Vu(r,i)||). (55) 

Here e(|| Vu(r, i)||) is appropriately chosen so that it is edge sensitive and is essentially zero away from a numerical 
shock or "discontinuity" . This leads to another shock capturing equation 

+ V • F(^.) - V • [di(||Vu(r, t)\\)S/u{r,t)] 

+ e(||Vu(r,i)|i). (56) 

Appropriate choice of the shock enhancing functional will result in a stable numerical algorithm. Obviously, the von 
Neumann and Richtmyer 0| artificial viscosity term, Eq. (^ , is a special case of the present formulation. 
The diffusion equation can be derived from Fick's law for mass flux, 

hir,t) = -DiVuir,t) (57) 

with Di being a constant. From the point of view of kinetic theory, this is an approximation to a quasi homogeneous 
system which is near equilibrium. A better approximation can be expressed as a super flux 

jq(r,i) = -^i?,VV29w(r,t), ((7 = 1,2,...), (58) 

9 

where Dq are constants and higher order terms {q > 1) describe corrections to mass flux by the influence of inhomo- 
geneity in density distribution and of flux-flux correlations. The mass conservation leads to 

du(r,t) „ . , , , , 
^^ = -V-jq(r,<) + s(r,i) 

= ^V-[i?,VV2««(r,i)]-t-s(r,t), (g-1,2,-.-), (59) 



where s is a source term which can be a nonlinear function describing chemical reactions. Equation ( |59| ) is a generalized 
reaction-diffusion equation which includes not only the usual diffusion and production terms, but also super diffusion 
terms. A truncation at the second order super flux, j2(r, i) = —Di'Vu{r,t) — D2W'^u{r,t), leads to the expression that 
is important in many phenomenological theories, such as the Cahn-Hilliard equation and the Kuramoto-Sivashinsky 
equation. The latter have been used for the description of a number of physical phenomena, such as the Taylor- 
Couette flow, parametric waves and pattern formation in alloys, glasses, polymers, combustion and biological systems 

il- 

In a shock wave simulation process, the distribution of the physical quantity under study can be highly inhomoge- 
neous and/or oscillatory. Hence, the present shock capturing algorithm can be made more efficient by incorporating 
a shock (edge) sensitive super diffusion operator 

+ V • Fin) = ^ V • [d,{\\Vu{r,t)m^''u{r,t)] 

9 

+ e(||Vw(r,t)||), (g-1,2,...). (60) 
Here dq{u, |Vw|) are shock (edge) sensitive diffusion fimctions. In most applications, a truncation at (7 = 2 is sufficient 

+ V • F(^.) = V • [di(||Vu(r, t)\\)Vu{r,t)] 

+ V-[d2(||Vu(r,i)l|)VV2w(r,t)] 

+ e(||Vu(r,i)||). (61) 
This equation can be simplified further for practical implementation 
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^iA_L2 + V . F{u) = Tidl Vu(r, t)\\)VMr, t) 
+ T2{\\Vu{r,t)\\)V\ir,t) 

+ e(||Vu(r,i)||), (62) 

where Fi and r2 are to be appropriately chosen to capture the shock. UsuaUy, Fi is a positive function, while r2 is 
negative. 

From the point view of image processing, both operators and are high pass filters. In the Fourier represen- 
tation, operators and are proportional to oj'^ and oj'^ of the frequency uj. Hence, is more sensitive to low 
frequency oscillations while has a large filter response for high frequency oscillations. Although both and 
become band pass filters on a discrete grid, their filter responses have similar properties as in their abstract operator 
forms. 

Note that many of these anisotropic terms presented in this section resemble those expressions derived in the 
last section for the artificial pressure tensor. Part- Although the starting points for these two approaches are entirely 
different, the results are obviously connected, notation needs to be simplified and possible confusions are to be avoided. 
For these reasons, we use the acronym ADOR for both approaches. 



IV. NUMERICAL EXPERIMENTS 



In this section, the numerical schemes developed from the analysis of kinetic theory and image processing analysis 
are utilized for numerical computations involving shock waves. Parameter selections for the ADOR method are 
specified. Obviously, there are infinitely many ways to construct anisotropic diffusion functionals. For simplicity, we 
test the following choices of the ADOR coefficients: 

di = dl ln[(Vu)2 + 1], d2 = d] \n[{Vu)'^ + 1], e^e^{V- (63) 

for Eq. (|§), 

Ti = 7i(l^^.|)^ r2 = 71(1^^1)^ e = (64) 

for Eq. del), 

Fi=72ln[|u,| + 1], F2=7|ln[|M.| + l], e = 0, (65) 

for Eq. (|6|), 

Fi = 7f ln[ul + 1], F2 - 7l + 1], e - 0, (66) 

for Eq. (H), and 

Fi=74ln[||Vu|| + l], F2=72ln[||V^/|| + l], e = 0, (67) 
for Eq. (|6^). In fact, coefficients in x- and y-direction can be chosen separately 

=-'.'{ 5:1: II 

Note that these expressions have one feature in common, i.e., they are all of low order in the gradient of u. This 
feature allows a sharp change in the solution to be handled by a coarse grid. Otherwise, if there are higher-order 
gradient terms, such as the von Neumann and Richtmyer form the computational grid near the shock front has 
to be refined to reduce the flux amplitude. A finer grid, in turn, leads to another stability problem as dictated by the 
CFL condition. 

In the rest of this section, the performance of the abovementioned prescriptions is examined for a detailed choice of 
ADOR parameters d\,dl,e^ and 7] (z = 1,2,3; j = 1,2). A few standard problems are employed to test the validity, 
and to demonstrate the robustness of the present approach. These problems include Burgers' equation in one and two 
space dimension, the incompressible Navier-Stokes equation and Euler equation with periodic boundary conditions. 

For spatial discretization, we use the discrete singular convolution (DSC) algorithm which was proposed as a 
potential approach for computer realization of singular convolutions p8|-^ . Mathematical foundation of the algorithm 
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is the theory of distributions [|^. Sequence of approximations to the singular kernels of Hilbert type, Abel type and 
delta type were constructed. Applications are discussed to analytical signal processing, Radon transform and surface 
interpolation. Numerical solutions to differential equations are formulated via singular kernels of delta type. By 
appropriately choosing the DSC kernels, the DSC approach exhibits global methods' accuracy for integration and local 
methods' flexibility for handling complex geometries and boundary conditions. Unified features of the DSC approach 
to differential equations were analyzed in details. In particular, we demonstrated |4l|] that different implementations 
of the DSC algorithm, such as global, local, Galerkin, collocation, and finite difference, can be deduced from a single 
starting point. The DSC algorithm is validated for the numerical solution of the Fokker-Planck equation [Q, the 
Schrodinger equation and the Navier-Stokes equation [|5|,|6j. It was also utilized to integrate the (nonlinear) 
sine-Gordon equation with the initial values close to a homoclinic orbit singularity ]39|] , for which conventional local 
methods encounter great difficulties and numerically induced chaos was reported for such an integration | |57[ |. The 
reader is referred to |Q for a detailed discussion of the method. The DSC parameters used in the present work are 
ct/A — 3.2 and M = 31 for the DSC kernel of regularized Shannon [Q. For time discretization, the explicit 4th-order 
Runge-Kutta scheme is used for Burgers' equation. The Navier-Stokes equation is integrated by using the implicit 
Euler scheme in association with a discrete singular convolution-alternating direction implicit (DSC-ADI) algorithm 
|p8[ . Further information is given in subsections below. 



A. Burgers' equation in one dimension 

Burgers' equation |^ is an important model for the understanding of physical flows. It appears customary to test 
new schemes in computational fluid dynamics by applying them to Burgers' equation. Despite much of the effort, 
numerical solution of Burgers' equation is still not a trivial task, particularly at very high Reynolds numbers where the 
nonlinear advection leads to shock waves. In fact, many standard computational algorithms fail to predict Burgers' 
inviscid shocks. 

Burgers' equation is given by 

du du 1 d'^u 

at+"9^"R^9^' ^ ' 

where u{x,t) is the dependent variable resembling the flow velocity and Re is the Reynolds number characterizing the 
size of the viscosity. The competition between the nonlinear advection and the viscous diffusion is controlled by the 
value of Re in Burgers' equation, and thus determines the behavior of the solution. We consider Eq. (|9|) using the 
following initial and boundary conditions 

u(x, 0) = sin(7ra:), 

u{0,t) = u{l,t) = 0. (70) 

Cole has provided an exact solution |Q for this problem in terms of a series expansion which is readily computable 
roughly for the parameter Re < 100. For the parameter Re = 100, the present calculations use 41 grid points in the 
interval [0,1] with a time increment of 0.01. Both Li and L^o errors at 11 different times are listed in TABLE I. 

We next consider inviscid Burgers' equation with same initial and boundary condition as given in Eq. (|^. In this 
case, the solution quickly develops into a sharp shock front a.t x — 1. The DSC algorithm does not work along because 
severe oscillations eventually turn into an uncontrolled error growth. This problem is treated by using the anisotropic 
diffusion approach discussed in the previous sections. The performance of four different prescriptions given in Eqs. 
(|63|)-(|66|) is examined. ADOR parameter selections are compared in FIG. 1 at 4 different times (t=0. 3, 0.5, 0.8, 2.0). 
These results are all obtained by using 101 grid points with a time increment of 0.002. 

In FIG. le, results for all the four different forms of coefficients are plotted for a detailed comparison. For the region 
away from the shock front, there is essentially no difference among four different forms. Although, results computed 
from different forms are slightly different at the shock front, four different choices of the diffusion coefficients have 
similar behavior and all are capable of correctly simulating Burgers' shock. 



The super diffusion coefficient chosen in Eq. (65) does not work very well as the shock front is distorted (see FIG. 
If). However, we have found that a combination of anisotropic diffusion and super diffusion does work better than the 
choice of a single term (results not shown). Similar results were also obtained by using 50 grid points. It is mentioned 
that the tests on the use of the edge enhancing functional does not result in a stable solution with the mesh system 
and time increment, chosen the present study. 
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B. Burgers' equation in two dimensions 



Let us consider Burgers' equation of the form 



Ut + UUx + VUy = -^{Uxx + Uyy) (7l) 



Re 
1 

Re' 



Vt + UV^ + VVy = —{Vjcx + Vyy) (72) 



in a square [0, 1] x [0, 1] with the initial values 



u{x, y, 0) — sin(7rx) sin(7ry) (73) 
y, 0) = [sin(7rx) + sin(27ra;)] + [sin(7r?/) + sin(27r?/)] (74) 



and boundary conditions 



u(0, y, t) = u(l, y, t) = u{x, 0, t) = u{x, 1, i) = (75) 

v(0, y, t) = v(l, y, t) = «(x, 0, t) = «(x, 1, t) = 0. (76) 

This problem has no analytical solution and is chosen to demonstrate that the present algorithm can be efficient, even 

with a very coarse mesh. This case is computed by using an ADOR parameter of 75* = 0.0006 with 41^ points. The 
contours of velocity field components at t=l are plotted in FIG. 2. 

C. The incompressible Navier-Stokes equations 

To test the present approach for shock capturing further we consider the Navier-Stokes equation 

Ut + UUx + VUy = + -^{UXX + Uyy) (77) 

Vt + UVx + VVy ^ -Py + -^{Vxx + Vyy) (78) 



with equation of continuity 



Re 



0, (79) 



where {u, v) is the velocity vector, p is the pressure. Re (Re> 0) is the Reynolds number and Re= oo defines the Euler 
equation. The domain of problem is a square [0, 2tt] x [0, 27r] with periodic boundary conditions. With appropriate 
initial values, the Euler equation can be used to describe a flow field of vertically perturbed horizontal shear layers 
around a jet. Bell et al ||] studied this case by a second order projection method. Recently E and Shu have 
employed this example to demonstrate the success of their high order ENO scheme for resolving the fine vorticity 
structure of the double shear layers. 

Case 1: Analytically solvable initial values 

The Navier-Stokes equation is analytically solvable for appropriate initial values 

u{x, y, 0) = — cos(a:) sin(y) 

v{x, y, 0) = sin(a::) cos(y). (80) 

The exact solution for this case is given by 

u{x, y,t) — — cos(x) sin(y)e~ 
v{x, y, t) ~ sin(a;) cos(y)e~ 

p(x,y,t) = -i[cos(2a;) +cos(2y)]e-*. (81) 

This provides a benchmark test for potential numerical methods in fluid dynamics. The implicit Euler scheme is used 
for the time integration and the DSC algorithm is utilized for the spatial discretization. The accuracy and reliability 



of this combination was previously tested |55| for this problem using a standard LU decomposition algorithm for 
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solving linear algebraic equations. Here, DSC-ADl algorithm js^l is used. We choose a grid of 32^ for the present 
calculation with a time increment of 0.001. The DSC-ADI results are summarized in TABLE II. Note that, for the 
inviscid case (Re=oo) the present result is accurate to the machine precision. It is evident that the accuracy of the 
DSC approach is extremely high, particularly when the Reynolds numbers are very large. 
Case 2: The Euler equation 

We now test our anisotropic diffusion approach for the Eulcr equation (Re= oo) with sharply varying initial values. 
This example is chosen to illustrate the ability of the present approach, for providing very fine resolution with a 
relatively coarse grid. The initial values are that of a jet in a doubly periodic geometry 

f tanh ( ^] 
\tanh(-^ 

w(a;,y,0) = (5sin(a:), (82) 

where 5 = 0.05 is used for the convenience of comparison with the previous study ||l^. This initial value describes 
the flow field consisting of horizontal shear layers of finite thickness, perturbed by a small amplitude vertical velocity, 
making up the boundaries of the jet. However, this problem is not analytically solvable. A pioneer work in this study 
was given by Bell et al in which they utilized a second-order Godunov scheme in association with a projection 
approach for divergence-free velocity fields with general boundary condition. With a periodic boundary condition, E 
and Shu have shown that their high-order ENO scheme ||l^ performs well for this problem. 

We first consider the parameter p = 7r/15, a case studied by Bell et al Q using a projection method with three 
sets of grids (128^, 256^ and 512^). E and Shu ||l^ computed this case by using both spectral collocation code with 
512^ points and their high order ENO scheme with 64^ and 128^ points. The spectral collocation code produced an 
oscillatory solution at i = 10 (see FIG. 1 of Ref. jl^), while the high order ENO scheme produced a defect at i = 6 as 
the channels connecting the vorticity centers are slightly distorted (see FIG. 2 of Ref. |l^). In the present simulation, 
we choose a 64^ grid for the computational domain with a time increment of 0.002. The prescription in Eq. ( |66| ) is 
used with the ADOR parameter of 7^ — 0.0006. The results at different times (t=4,6,8,10) are plotted in FIG 3. It 
is seen that our solutions are smooth (some non-smooth features in the contour plot is due to the fact that the grid 
is very coarse) and stable for this case. In particular, no distortion is found in vorticity contours at t=6. For early 
times, present results compare extremely well with those of the spectral collocation code computed with 512^ points. 
There are no spurious numerical oscillations during the entire process. 

Finally, we perform simulations by the present approach using the discontinuous initial data {p 0) as in . The 
evolution under the Euler equation leads to a periodic array of large vortices, with the shear layers between the rolls 
being thinned by the large straining field. It is known that, for the present choice of parameters the solution quickly 
develops into a roll-up process with smaller and smaller scales, and the resolution is lost eventually with a fixed grid 
for local methods. We compute this case by using 128^ grid points as in jl^. A time increment of 0.002 is used. The 
ADOR parameter is chosen as 75* = 0.0006. The present results at different times (t=4,6,8,10,12,14) are plotted in 
FIG 4. 



V. CONCLUSIONS 



Connection is made between digital image processing and computational fluid dynamics. The evolution of an 
image surface under a partial differential operator can be viewed as a form of image processing. Computationally, 
numerical shock capturing can be formulated on the lines of iterative edge-detection. Hence, techniques developed in 
the computational fluid dynamics can be used for image processing and vice- versa. This paper introduces the method 
of anisotropic diffusion oscillation reduction (ADOR) , an approach which has its roots in image processing, for shock 
wave computations. 

In fact, the ADOR method is much similar to the artificial viscosity algorithm. Physical origins and mathematical 
properties of the artificial viscosity are discussed from the kinetic theory point of view. The form of pressure tensor is 
derived from the first principles of quantum mechanics. Quantum kinetic theory is utilized to arrive at macroscopic 
transport equations from the microscopic quantum theory. Macroscopic symmetry is used to simplify the phenomeno- 
logical pressure tensor expressions. The latter provides the basis for the design of artificial viscosity. The original von 
Neumann-Richtmyer form of artificial viscosity fits into a special case of the present generalizations. 

The anisotropic diffusion, which is essentially an image processing technique, is modified for shock capturing. The 
technique preserves image edges by introducing little diffusion, where the image gradient is large, while it provides 
a substantial diffusion coefficient at smooth parts of the image. In the present shock capturing algorithm, the edge- 
sensitive diffusion is introduced at a shock front so that large oscillations can be efficiently eliminated. An edge 
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enhancing functional and edge-detected super diffusion operators are proposed for shock capturing. These terms are 
introduced from the general structure of conservation laws. In fact, kinetic theory analysis and anisotropic diffusion 
argument lead to a number of similar expressions for shock wave treatments. Hence, the acronym ADOR is referred 
for both approaches. 

The reliability and robustness of the ADOR method is explored in association with the discrete singular convolution 
(DSC) algorithm [^-^. The DSC algorithm was previously validated by handling many linear and nonlinear 
problems and is a potential algorithm for the computer realization of some singular convolutions. A few detailed 
prescriptions for the anisotropic diffusion functional are considered in this work. The coefficients in these prescriptions 
are chosen to be of lower order in gradient so that a coarse grid can be used. A number of standard test examples, 
including (inviscid) Burgers' equation in one and two spatial dimensions, the incompressible Navier-Stokes and the 
Eulcr equations, are employed for the present test computations. 

For Burgers' equation, we have examined the accuracy for spatial and temporal discretizations. A few ADOR 
coefficients given by Eqs. (|63|)-(|66|) are tested. The anisotropic diffusion term is found to be very robust in association 
with many edge-detected coefficients, while the super diffusion operator produces an over shot at Burgers' shock front. 
The edge enhancing functional apparently does not perform well for this problem. The ADOR approach is found to 
be very stable for the 2D inviscid Burgers' equation even with a very coarse grid. 

For the incompressible Navier-Stokes equation, the ADOR approach is used in association with a DSC-ADI algo- 
rithm. The accuracy of the algorithm was tested by an analytically solvable case and the machine precision is found 
at the inviscid limit of the Navier-Stokes equation. The algorithm was then applied to the study of doubly periodic 
shear layers, which was chosen to demonstrate the capability of the ADOR method for more difficult problems. These 
results indicate that the proposed method has the potential to capture shock waves. Work on the implementation of 
the present approach to more complex fluid flow systems, including general boundary conditions and for compressible 
flow is under progress. 
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TABLE I. Li and Lao errors of the DSC solutions for Burgers' equation 



Time 


Li 




0.4 


2.1(-4) 


2.8(-3) 


0.8 


2.8(-4) 


4.4(-3) 


1.2 


4.2(-5) 


6.8(-4) 


1.6 


5.8(-6) 


9.4(-5) 


2.0 


9.4(-7) 


1.2(-5) 


3.0 


3.8(-8) 


4.0(-7) 


5.0 


2.1(-10) 


2.2(-9) 


10 


1.5(-11) 


4.0(-ll) 


20 


3.2(-12) 


5.3(-12) 


30 


1.1(-12) 


1.8(-12) 


40 


4.1(-13) 


6.4(-13) 



TABLE II. L2 errors of the DSC-ADI solutions for the 2D Navier-Stokes equation computed with 32 grid point in each 
dimension 



Re 


t=l 


t=2 


t=3 


t=4 


u 


V 


u 


V 


u 


V 


It 


V 


20 


6.1(-11) 


6.1(-11) 


l.l(-lO)'' 


l.l(-lO) 


1.5(-10) 


1.5(-10) 


1.8(-10) 


1.8(-10) 


10^ 


1.2(-12) 


1.2(-12) 


2.3(-12) 


2.3(-12) 


3.4(-12) 


3.4(-12) 


4.4(-12) 


4.4(-12) 


10^ 


9.5(-13) 


9.3(-13) 


1.9(-12) 


1.8(-12) 


2.8(-12) 


2.8(-12) 


3.8(-12) 


3.7(-12) 


10"* 


8.4(-13) 


8.3(-13) 


1.7(-12) 


1.7(-12) 


2.6(-12) 


2.6(-12) 


3.5(-12) 


3.5(-12) 


10^ 


4.6(-13) 


4.6(-13) 


9.7(-13) 


9.7(-13) 


1.4(-12) 


1.4(-12) 


1.9(-12) 


1.9(-12) 


10^ 


6.3(-13) 


6.3(-13) 


1.2(-12) 


1.2(-12) 


1.7(-12) 


1.9(-12) 


2.6(-12) 


2.6(-12) 


00 


1.5(-15) 


1.5(-15) 


1.6(-15) 


1.6(-15) 


1.6(-15) 


1.6(-15) 


1.7(-15) 


1.7(-15) 



L2 = 9.1(-04) and ^ L2 = 4.9(-04) were reported in Ref. |T||. 
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Figure Captions 

FIG. 1. A comparison of various ADOR cocfiicicnts for tlie numerical solution of inviscid Burgers' equation, (a) 
d\ = 0.0009, d\ = 0, el = 0; (b) 7! = 0.0022,7] = 0; (c) 7? = 0.0018, 7I = 0; (d) 7? = 0.0015, 7I = 0; (e) Comparison 
of (a), (b),(c) and (d); (f) 7? = 0,7| = -98 x 10"^ 

FIG. 2. The ADOR solution for the 2D inviscid Burgers' equation with 41^ points, 7^ = 0.02, 7I = 0. Left: the u 
field; right: the v field. 

FIG. 3. The vorticity contours for the 2D Eulcr equation by the ADOR method with 64^ points, 7^ = 0.0006. Up 
left: t=4; up right: t=6; low left: t=8; low right: t=10. 

FIG. 4. The vorticity contour for the 2D Euler equation, with discontinuous initial data, by the ADOR method 
with 128^ points, 7^ = 0.0006. Up left: t=4; up right: t=6; middle left: t=8; middle right: t=10; low left: t=12; low 
right: t=14. 
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FIG. 1. A comparison of various ADOR coefficients for the munerical solution of inviscid Burgers' equation, (a) 
dl = 0.0009,4 = 0,ei = 0; (b) = 0.0022, = 0; (c) 7^ = 0.0018, 7I = 0; (d) 7? = 0.0015, 7I = 0; (e) Comparison 
of (a), (b),(c) and (d); (f) = 0,7| = -98 x 10"^ 
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FIG. 2. The ADOR solution for the 2D inviscid Burgers' equation with 41^ points, jf = 0.02, = 0. Left: the u 
field; right: the v field. 
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FIG. 3. The vorticity contours for the 2D Euler equation by the ADOR method with 64^ points, 7^ = 0.0006. Up 
left: t=4; up right: t=6; low left: t=8; low right: t=10. 
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FIG. 4. The vorticity contour for the 2D Euler equation, with discontinuous initial data, by the ADOR method with 
128^ points, 7f = 0.0006. Up left: t=4; up right: t=6; middle left: t=8; middle right: t=10; low left: t=12; low right: 
t=14. 
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